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Abstract. Western European Water Frog species can be morphologically very similar, especially Pelophylax ridibundus 
and P. perezi. As these species, together with P. lessonae, produce hybridogenetic lineages with highly variable phenotypes, 
the morphological identification of taxa in this genus remains challenging. This is especially problematic in the current 
context of the spread of the invasive P. ridibundus complex in native populations across Western Europe. In this study, we 
present a novel molecular method for the identification of western European Water Frogs (genus Pelophylax) based on the 
PCR-RFLP genotyping of nuclear species-specific loci. Using a combination of restriction digestions of species-specific 
genetic markers, this method facilitates reliable specific identification in the Pelophylax perezi-grafi system (PG system) 
and the Pelophylax lessonae-esculentus system (LE system). Our method provides several advantages compared to other 
genetic techniques previously used. It relies on PCR amplification allowing non-invasive sampling followed by RFLP to 
assess diagnostic SNPs in the nuclear genome, rendering obsolete the need for sequencing. Lastly, it is easier to replicate 
than methods based on microsatellites. Importantly, this novel protocol allows the identification of all five taxa of Water 
Frogs that occur in Western Europe. 
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Introduction 


The western Palearctic Water Frogs (genus Pelophylax) 
constitute one of the best-known examples of hybridoge- 
netic species complexes, giving rise to several hybridoge- 
netic systems across Europe (GUNTHER 1990). These com- 
plexes are composed of a fixed bisexual hybrid, originating 
from hybridization events between Pelophylax ridibundus 
and another species, usually coexisting with at least one of 
its parental species. The hybrid taxa are maintained by a 
hemiclonal reproductive system where, during gametogen- 
esis, the hybrid taxon systematically excludes the genome 


of one of the parental species, frequently the genome from 
the parent species that lives in sympatry with it. Pre-meiotic 
DNA synthesis is followed by a compensatory duplication 
of the remaining genetic material, which generates gametes 
that contain the genome of just one of the parental species 
(VINOGRADOV & CHUBINISHVILI 1999). Thus, recombina- 
tion usually does not take place except in triploids (CHrRs- 
TIANSEN & REYER 2009) and the genome that is included in 
the gametes of the hybrid taxon (typically but not always a 
P. ridibundus hemigenome) is transmitted clonally. Conse- 
quently, only the somatic cells of the hybrid contain genetic 
material of both parents, while gametes are clones of a sin- 
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gle parent’s genetic material. Given this particular repro- 
ductive system, hybrid maintenance relies on backcrossing 
to the parental species whose genome is eliminated: cross- 
ings between the hybrid taxa and the co-occurring parental 
species produce hybrid taxa, while crossings between in- 
dividuals of the hybrid taxa will typically fail, as the clon- 
ally transmitted genome accumulates too many deleteri- 
ous mutations to function in a diploid state (GRAF & POLLS 
PELAZ 1989, DUFRESNES & MAZEPA 2020). The hybrid taxa 
are thus the avatars of hemiclones that reproduce by “para- 
sitizing” a sexual species; these hemiclones function well 
when combined with the sexually reproducing genome of 
another species, but are usually unable to persist on their 
own (see DUFRESNES & MAZEPA 2020 for a review). 

The most common and most studied system is the 
Pelophylax lessonae-esculentus system (the LE system, 
BERGER 1973, HEPPICH 1978), found throughout Europe 
from western France to the Volga-Urals region of Russia 
(GRAF & POLLS PELAZ 1989, Gasc et al. 2004). This spe- 
cies complex, composed of P. kl. esculentus, a fixed bisexual 
hybrid resulting from hybridization events between P. ridi- 
bundus and P. lessonae, and the parental species P. lessonae, 
upon which the hybrid taxon depends for its reproduction 
(BERGER 1973, HEPPICH 1978, GRAF & MUELLER 1979, UZ- 
ZELL et al. 1980), has allowed the characterization of the 
particular reproductive system of the Water Frog complex- 
es. Several other hybridogenetic systems involving these 
three taxa are distributed across Europe (GRAF & POLLS 
PELAZ 1989, DUFRESNES & MAZEPA 2020). They differ in 
population composition, especially with regard to the spe- 
cies coexisting with the hybrid, and in the genome trans- 
mitted by the hybridogenetic taxon (ridibundus in most 
systems, including the LE system, but rarely lessonae in the 
ridibundus-esculentus RE system, see details in DUFRESNES 
& MAZEPA 2020). They also differ in the ploidy of the hy- 
brid P. kl. esculentus (mainly diploid individuals with one 
P lessonae and one P. ridibundus genome in the LE system 
or including triploids with two genomes of either P. lesso- 
nae or P. ridibundus in other systems) and of the gametes 
produced by esculentus (mainly haploid in the LE system, 
often diploid in some other systems). Lastly, some pure 
populations of the P kl. esculentus hybrid taxon are known 
to persist through the coexistence of LR diploids with LLR 
or LRR triploids (the esculentus-esculentus EE system). In 
some of these systems recombination and sexual repro- 
duction have even been restored, giving rise to sexually re- 
producing populations of the hybrid taxon P kl. esculentus 
(see CHRISTIANSEN & REYER 2009). 

Another taxon of hybrid origin is Pelophylax kl. grafi, a 
fixed bisexual hybrid between Pelophylax perezi and P. ridi- 
bundus. Coexisting with P. perezi to form the perezi-grafi 
system (PG system), it is commonly encountered in the 
Mediterranean and Atlantic regions of southern and south- 
western France and in northeastern Spain (GRAF et al. 1977, 
UZZELL & TUNNER 1983, GRAF & POLLS-PELAZ 1989, HOTZ 
et al. 1994, CROCHET et al. 1995, PAGANO et al. 20014, b, DAF 
et al. 2006). The current consensus is that natural hybrido- 
genetic populations were originally constituted exclusively 


of individuals of P. perezi and P. kl. grafi, except in the re- 
stricted area of overlap with the LE system (see Dar et al. 
2006). The current widespread occurrence of P. ridibundus 
within the distribution of the PG system in France is as- 
sumed to be the result of expansions from spots of recent 
human-mediated translocations that triggered the inva- 
sion of this region by P. ridibundus (PAGANO et al. 2001, b, 
2003, HOLSBEEK & JOORIS 2010, HOLSBEEK et al. 2010, Du- 
FRESNES et al. 2017a). As in the LE system, the hybrid taxon 
is thought to persist by backcrossing with P perezi and by 
transmitting only the P ridibundus genome. 

The origin of P. kl. grafi remains speculative and two 
main hypotheses are still under discussion (ARANO et al. 
1995). One hypothesis proposes direct hybridization be- 
tween P perezi and P. ridibundus during past periods of 
range overlap. Another hypothesis involves hybridization 
between P perezi and P. kl. esculentus (with the latter acting 
as the donor of the P. ridibundus hemigenome) where the 
distributions of the LE and PG systems currently overlap 
in western France (PAGANO et al. 20012, b, DAF et al. 2006). 
Alternatively, P. ridibundus could be native (and not intro- 
duced as is widely believed) to southeastern France, where 
it is currently the only Pelophylax in most of Provence and 
Alps regions (LESCURE & DE MASSARY 2012, own unpub- 
lished data) and comes into contact with the PG system 
along the lower Rhône valley. However, Horz et al. (1994) 
argued that the low clonal diversity of P. kl. grafi compared 
with P kl. esculentus and the fact that the system was re- 
stricted to the eastern limits of the geographic distribution 
of P. perezi supported the idea that P. kl. esculentus is the 
source of the P. ridibundus genome in the PG system. 

One of the major issues precluding a better understand- 
ing of the ecology and evolution of the European Wa- 
ter Frogs is the difficulty of reliable species identification 
(GUNTHER 1990, PAGANO & JOLY 1999). In the systems in- 
volving F kl. esculentus, morphology and bioacoustics dif- 
fer between the parental species and the hybrid (GUNTHER 
1990, PLOTNER 2005), but gene dosage effects in triploid 
hybrids that contain a double genetic makeup, either from 
P. ridibundus or from P. lessonae, blur the validity of mor- 
phological and acoustic delimitations (PLOTNER 2010, 
HOFFMANN & REYER 2013, MAYER et al. 2013). In the case 
of P. kl. grafi, phenotypic differentiation is even more dif- 
ficult, since P. perezi and P. ridibundus are extremely simi- 
lar (GUNTHER 1990) and bioacoustic differences have not 
been well established yet. Therefore, we still do not know 
how to confidently identify P. kl. grafi based on morpho- 
logical or acoustic characters. 

Accurate identification of interspecific hybrid taxa and 
parental species is not only crucial for understanding the 
complex reproductive system of hybridogenetic species, 
but it also represents an important issue in conservation 
biology (e.g., CONGIU et al. 2001). Several studies suggest 
that the presence of P ridibundus in western European 
countries (e.g., in the south of France, Spain, Switzerland) 
is the result of multiple human-mediated translocations 
(ARANO et al. 1995, PAGANO et al. 1997, PAGANO et al. 20014, 
2003, VORBURGER & REYER 2003, ZEISSET & BEEBEE 2003, 
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HOLSBEEK & JOORIS 2010, HOLSBEEK et al. 2008, 2010). 
Pelophylax ridibundus is currently responsible for the de- 
crease of native Water Frog species in many areas and it is 
considered a threat to other amphibian species in the re- 
gion (ARANO et al. 1995, PAGANO et al. 2003, VORBURGER 
& REYER 2003, SCHMELLER et al. 2007, HOLSBEEK & JOORIS 
2010, HAUSWALDT et al. 2012, LEUENBERGER et al. 2014). 
Therefore, molecular identification methods of Water Frog 
complexes are essential tools for monitoring the spread of 
the invasive P. ridibundus, for establishing conservation 
measures, as well as for ecological and evolutionary studies. 

In European Water Frogs, interspecific mitochondrial 
gene transfer (SPOLSKY & UZZELL 1984, 1986, PLOTNER et 
al. 2008) and the presence of hybridogenetic taxa make 
mitochondrial markers unreliable, because they are mater- 
nally inherited and thus cannot identify hybrid taxa. Such 
problems do not occur with nuclear markers, and the mo- 
lecular identification of Water Frogs has long been based 
on allozymes (UZZELL & BERGER 1975, UZZELL & Hotz 
1979, HOTZ & UZZELL 1982, GUNTHER & PLOTNER 1994, 
Horz et al. 1994, BUCKLEY et al. 1994, BEERLI et al. 1996). 
However, this method comes with technical challenges 
such as the need for large amounts of tissue and the ne- 
cessity to preserve samples at very low temperatures upon 
collection; it also suffers from poor repeatability, as al- 
leles are defined by their relative migration speed, render- 
ing identification of specimens without reference samples 
difficult. Nuclear DNA-based methods that alleviate some 
of these difficulties have been proposed, using microsatel- 
lites (GARNER et al. 2000, ZEISSET et al. 2000, Hotz et al. 
2001, CHRISTIANSEN 2005, SANCHEZ-MONTES et al. 2016, 
SaGonas et al. 2020), Southern blot analysis of enzyme- 
digested genomic DNA (ToGNaARELLI et al. 2014), or ran- 
dom amplified polymorphic DNA markers (RAPD, ZEIs- 
SET & BEEBEE 1998) as alternatives. Microsatellites, how- 
ever, suffer from a lack of diagnostic loci and usually re- 
quire assignment analyses in large samples to allow species 
identification. RAPD approaches suffer from low repeat- 
ability, while Southern blot methods, which do not rely on 
DNA amplification, require large amounts of tissue. None 
of these methods is thus suitable for non-invasive specific 
identification of single specimens or are easily transferable 
from one lab to another. 

Alternative methods based on PCR-amplification of 
diagnostic nuclear loci have thus been proposed. A PCR- 
RFLP method by Parre te et al. (2011) and a PCR-mi- 
gration method by Hauswa pr et al. (2012) have been 
developed to identify the Water Frogs involved in the LE 
complex (P ridibundus, P. lessonae and P. kl. esculentus), 
with the latter method having been successfully adapted 
to separate P. ridibundus from P. shqipericus (Vucic¢ et al. 
2018). However, these methods have not been tested in 
the Pelophylax perezi-grafi (PG) system and are thus for 
the moment unsuitable in large areas of Western Europe. 
For the PG system, the only DNA-based method current- 
ly available uses a set of 16 microsatellite loci (SÁNCHEZ- 
Montes et al. 2016); like allozymes, microsatellites can 
suffer from poor transferability between labs, and identi- 
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fication of specimens without reference samples is diffi- 
cult. Moreover, given that in some regions the LE system 
and the PG system overlap, a single reliable method for the 
identification of the three parental species and the two hy- 
brid taxa involved in the PG and LE systems is necessary. A 
technique based on nuclear markers with biparental inher- 
itance and where genotyping is possible without the need 
for reference samples will provide a crucial tool for the reli- 
able identification of the different Water Frog taxa in these 
two species complexes. 

In this study we present a fast, easy and non-invasive 
PCR-RFLP method that allows the separation of the five 
widespread European Water Frog species: Pelophylax ridi- 
bundus, P. perezi, P. lessonae (including bergeri, which is 
sometimes treated as a valid species) and the hybridogenet- 
ic P kl. esculentus and P. kl. grafi. Because recombination 
in the hybridogenetic lineages cannot be totally excluded 
(CHRISTIANSEN & REYER 2009), we aimed at providing at 
least two independent markers to diagnose each taxon. We 
first selected nine nuclear loci found to be informative at 
interspecific level by previous studies or already published 
for Pelophylax species identification, including the two loci 
used by PATRELLE et al. (2011) and Hauswarnpr et al. (2012, 
Table 1). For those where amplification was successful, se- 
quences were generated for reference samples of the three 
parental species (P. ridibundus, P. perezi and P. lessonae, 
see Table 2). This allowed us to identify five nuclear genes 
with potentially diagnostic substitutions between at least 
two of the parental species (Table 3). After screening po- 
tential enzymes with restriction sites around the diagnos- 
tic positions, we tested enzymatic digestions on large series 
of samples representing the five target taxa to validate the 
diagnostic character of these substitutions. Because most 
of the samples used in this last step had not been identi- 
fied previously by other methods, we employed at least two 
genes to identify each sample and assessed the concord- 
ance between the identification provided by each locus 
to validate the identification. This concordance between 
markers also provides a test for the reliability of the diag- 
nostic sites and allows detecting rare recombination events 
(CHRISTIANSEN & REYER 2009). 


Material and methods 
Sampling 


We analyzed a total of 1,111 individuals belonging to all five 
Western European Water Frog species, Pelophylax ridibun- 
dus, P. lessonae, P. Kl. esculentus, P. perezi and P. kl. grafi 
(Supplementary File 1). Water Frog sampling in France had 
been authorized by the following permits from competent 
regional administrations: AP 02/2007 & 03/2007 by the 
Préfecture de la Vendée; AP 07-2549 by the Préfecture de 
la Sarthe; AP by the Préfecture des Ardennes in 2007; AP- 
2013-06; AP-2013 274-0002; AP-2013-325.0010; AP-2013- 
189.0001; AP 2019-s-07 by the Préfecture de Hérault et 
du Gard, and samples from other countries were retrieved 
from the BEV tissue collection. 
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Table 1. List of the nine nuclear genes evaluated in this study with their corresponding primers’ names, sequences and sources. An- 
nealing temperature is only given when amplification was successful. The “Result” column identify which genes were selected. For the 
selected genes, the last column reports the restriction enzymes that were retained after our tests. * Details for the digestion of NTF3 
gene by the Xap] enzyme are given in Supplementary Files 2 (PCR-RFLP method) and 3 (alignment). 


Annealing Restricti 
Gene Primers Sequence 5’ - 3’ Source tempera- Result esthietion 
enzyme 
ture 
CXCRA ; j 
haon CXCR4_Rana.F 5’-TTCACCCTTCCATTCTGGTC-3 NEWMAN et © Notselected(No | 
CXCR4_Rana.R 5'-GCCACGGCTTCTGTGATAG-3 al. (2016) amplification) 
receptor type 4) 
sia 5’-TCGAGTGCCCCGTGTGYTTYGAY- 
: TA-3’ Frost et al. Not selected(No 
Soleo SIAI (F)SIA2 (R) 5: Ga AGTGGAAGCCGAAGCAGSWYT- (2006) = amplification) 7 
GCATCAT-3’ 
irs? 5’-GGATCACTCGGCTCGTGCGTCGAT- Nöt selectediio 
(Int It ITS2-hinITS2-riick GAAG-3 PATRELLE et ification i 
ze eae UME SATUS 5°-CGGGGATTCGGCGCTGGGCTCTTC- al. (2011) — ee 
scribed spacer 2) CC?’ perezi or graft) 
SAI-1 5’- TCCATACAAATGTGCTAAGTAG- Not selected 
_  PEL-SA-F1PEL- > HAUSWALDT 3 (Length differences 
(Serum albumin SA-R2 Gri et al. (2012) ete not entirely diag- 
intron 1) 5’-GACGGTAAGGGGACATAATTCA-3’ i nae y Mag 
5’- TCTTCCTTATCTTTGTGGCATCCAC- Not selected" (not 
NTF3 GCTA-3’ NEWMAN et a : : 
(Neurotrophin 3) NTF3-FNTF3-R © aCaTTGRGAATTCCAGIGTTTGTC- al. (2016) 61°C fully diagnostic Xapl(RzAATTY) 
; site) 
GTCA-3 
RAGI 5’- AGCTGCAGYCAGTACCA- Dra II(RG- 
Peal al aaa RR CAAAATG-3° one st 59°C Selected zGNCCY)BSM I 
brotein 5 5’- AATTCAGCTGCATTTCCAATGTC-3’ ™ (GAATGCNz) 
TYR TYR1bRana (F) 5’-AGGTCCTCTTRAGCAAGGAATG-3’ NEWMAN et dC. Seeded Bal I(TGGzC- 
(Tyrosinase 1) TYR1gBufo(R) 5°-TGCTGGGCATCTCTCCAGTCCCA-3’ al. (2016) CA) 
pöüe POMC SH Rana Ž-ATAYGTCATGAGCCACTTCCGCTG- < Ysanne l 
(Pro-opiomelano- FIPOMC_SH. GAAS HAUSWALDT 61°C Selected Line OTIZ 
tin A) Rany R — Ž-GTCTTTGGGTGGWCCTTC- ) RAC) 
corun ana_. CCATCG-3’ pers. com. 
CMYC? DUBEY & 
(Cellular myelo- CMYC2 MYEL F1 5’-CAGTGAATGACAGCATTTCCAG-3’ : 
cytomatosis intron CMYC2 MYEL R3 5’-GTCAAAGCCTTCAAAGACCATTG-3’ aa ae SSelected aca 


2) 


DNA sampling and extraction 


Two types of samples were used in this study: 1) mostly 
unidentified tissue samples stored in 95% ethanol at -20°C 
in the tissue collection of the Biogéographie et Ecolo- 
gie des Vertébrés (BEV) team from the UMRs5175 - Cen- 
tre d’Ecologie Fonctionnelle et Evolutive (Montpellier, 
France); these samples include both toe clippings or buc- 
cal swabs from live frogs (not preserved after sampling) 
or tissue samples from vouchers previously preserved in 
95% ethanol in the BEV voucher collection; and 2) a se- 
ries of 83 samples from the five species originating from 
the Ardennes (northeastern France), southeastern and 
southwestern France, and northeastern Spain, collected 
between 1996 and 2008, stored in ethanol in the PAGA- 
No and PATRELLE tissue collection by the team “EA 4688 
USC ANSES - VECPAR’” (Université Reims Champagne- 
Ardenne, France). These 83 samples had been previous- 
ly identified using specific allozymic markers or a PCR- 


RFLP-based method on ITS2 (PAGANO et al. 1997, 20018, 
b, Dar et al. 2006, PATRELLE 2010, PATRELLE et al. 2011). 

In a first step, we selected a series of P. ridibundus, 
P. lessonae and P. perezi samples from the BEV sample col- 
lection, originating from areas where no other species are 
known to occur. These reference samples and samples from 
the PaGaNo and PaTRELLE collection, previously identified 
with other methods, were used to screen nuclear loci for 
candidate diagnostic single-nucleotide polymorphisms 
(SNPs). In a second step, a large series of DNA samples 
from the BEV and the PaGANo and PaTRELLE tissue collec- 
tions were used to assess the reliability of the method and 
to provide genetic identification to local naturalists and 
managers of protected areas (see details in Results and in 
Supplementary File 1). 

We used three different DNA extraction protocols. 1) A 
modification of the NaCl extraction protocol from Bru- 
FOLD et al. (1992), using an extraction buffer containing 
Tris (1M, pH 8), NaCl (5M), EDTA (0.5M, pH 8) and dH,O. 
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Samples were incubated in a mix of 410 ul of extraction 
buffer, 80 ul of SDS (10%), and 10 ul of Proteinase K at 65°C 
for 30 min before centrifugation at 13,000 rpm for 5 min. 
A volume of 180 pl of NaCl (5M) was added to the super- 
natant and centrifuged again at 13,000 rpm for 5 min. Fi- 
nally, 420 ul of cold isopropanol was added to the result- 
ing supernatant and a final centrifugation of 13,000 rpm 
for 5 min was performed. We washed the DNA pellet with 
250 ul of EtOH 80% before elution in 100 ul of water. 2) Sec- 
ondly, we used the Extract-N-Amp”™ Plant PCR Kit (Sigma- 
Aldrich, Saint-Louis, USA). Samples were mixed with 50 ul 
of E-Buffer, and after incubation at 95°C for 10 min, 50 ul of 
D-Buffer were added. 3) Lastly, some samples were extract- 
ed with the DNeasy Blood and Tissue Kit (Qiagen, Hilden, 
Germany) following the manufacturer’s protocol. All three 
extraction protocols showed good results for buccal swabs. 
However, for tissue samples, the NaCl extraction protocol 
was chosen given its better performance. 


Nuclear marker selection 


To identify loci allowing species identification, nine po- 
tential nuclear markers were selected based on previous 
studies, and primers for the amplification of these mark- 
ers were retrieved from the literature (see list in Table 1). 
Polymerase Chain Reaction (PCR) amplifications were op- 
timized by performing gradient PCRs in search of the opti- 
mal annealing temperature. PCR reactions were conducted 
in 20 uL volumes with 2 uL of DNA, 10 uL of 1X Taq Poly- 
merase (Extract-N-Amp™ PCR ReadyMix, Sigma-Aldrich, 
Saint Louis, USA), 0.5 uL of each primer (10 uM), and 7 uL 
of sterile water. The PCR amplification protocol was as fol- 
lows: a denaturation step of 5 minutes at 94°C; 41 cycles of 
denaturation at 94°C for 1 minute, annealing at specific an- 
nealing temperature (Table 1) for 1.5 minutes, elongation at 
72°C for 1 minute; a final elongation of 10 minutes at 72°C. 

Reference samples of the parental species were PCR-am- 
plified for each of the selected markers in Montpellier (BEV 
samples) or in Reims (PAGANO and PATRELLE collection) to 
identify potentially diagnostic SNPs (details in Table 2 and 
Supplementary File 1). Sequencing was done in both direc- 
tions with the PCR primers; samples from the BEV collec- 
tion or from the PAGANO and PaTRELLE collection were se- 
quenced by Eurofins Genomics (Ebersberg, Germany) or 
Genoscreen (Lille, France), respectively. Sequences were 
aligned using CodonCode Aligner V.4.2.3 sequence analy- 
sis software (CodonCode Corporation, Centerville, USA) 
or ClustalW (THompson et al. 1994) included in MEGA V.5 
(TAMURA et al. 2011). For the CMYC2 marker, we retrieved 
sequences from GenBank for several individuals of P ridi- 
bundus, P. perezi and P. lessonae (Table 2). Representative 
samples of newly generated sequences from each marker 
and parental species have been deposited in GenBank (see 
Table 2 and Supplementary File 1 for details). 

Sequences of each marker were inspected visually across 
all species using the alignment viewer in MEGA V.5 (Ta- 
MURA et al. 2011) to identify potentially diagnostic SNPs. 
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Table 2. Number of sequenced individuals for the four selected 
markers (or individuals retrieved from Genbank for CMYC2). 
Some individuals have been sequenced to check ambiguous re- 
striction patterns and have not been uploaded on GenBank. 


CMYC2 
RAG1 TYR1 POMC (GenBank) 
P. ridibundus 32 32 5 13 
P. perezi 14 9 
P. lessonae 6 6 6 8 


Among the six successfully sequenced markers (Table 1), 
five provided SNPs with alternative alleles in at least one of 
the three parental species: Recombination-activation pro- 
tein 1 (RAG1), Tyrosinase (TYR1), Pro-opiomelanocortin 
(POMC), Neurotrophin-3 (NTF3), and Cellular myelo- 
cytomatosis intron 2 (CMYC2). 

In addition, when the initial RFLP tests generated am- 
biguous patterns or discordance between the RAGi and 
TYRi markers for some samples, we sequenced some sam- 
ples for these markers to inspect the diagnostic SNPs di- 
rectly from the sequences. In each case, sequencing re- 
vealed the expected base for the species and removed the 
discordance between markers. After confirmation to that 
effect had been found in these first sequencing trials, sub- 
sequent conflicting or ambiguous results were resolved 
by repeating the RFLP steps until results were compatible 
across markers or the conflict between markers were con- 
firmed (see below). Most of these sequences were checked 
for their diagnostic positions, however, they were not cor- 
rected in their entirety, have not been deposited in Gen- 
Bank, and the respective individuals are identified as “se- 
quencing” in Supplementary File 1 instead of having Gen- 
Bank accession numbers. 


PCR-RELP based method: 
digestion by restriction enzymes 


Seeking to provide a quick and affordable method to iden- 
tify Water Frog specimens, we developed restriction frag- 
ment length polymorphism protocols (RFLP) to genotype 
samples at the candidate diagnostic SNPs identified in the 
previous step. We used the RestrictionMapper software 
(available from www.restrictionmapper.org) and the CLC 
sequence Viewer 6 software (Giagen, Aarhus, Denmark) 
to find restricting endonucleases (restriction enzymes) 
whose recognition sequence include the candidate SNPs. 
When several enzymes were available for the same SNP, we 
selected the ones that were cheaper and easier to obtain. 
We retained the following enzymes for further tests: DralI 
(= EcoO1o9I) and BsmI (= Mvai269I1) for digesting the 
RAG1 marker, Ball (= MscI) for the TYR1 marker, SspI for 
CMYC2, Xap] (= Apo I) for NTF3, and HinclI for POMC. 
All enzymes were ordered from the Thermo Scientific Fast- 
Digest Restriction Enzymes line (Thermo Fisher Scientific, 
Waltham, USA). 
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Table 3. Restriction patterns of the selected enzymes. For each marker (DNA fragment) and restriction enzyme, we provide the size 
of the fragment(s) resulting from the restriction digestion for each species. * The BsmI restriction site on RAG] is not fully diagnostic 
if Italian samples of the LE system are included (see Supplementary File 2 and 3). 


DNA fragment TYRI1 (646bp) 


RAGI (983bp) 


POMC (476bp) CMYC2 (298bp) 


Restriction enzyme Ball Drall 

P ridibundus 646 394/589 
P. perezi 305/341 983 

P. kl. grafi 646/305/341 983/394/589 
P. lessonae 646 983 

P kl. esculentus 646 983/394/589 


BsmI* Hmcll Sspl 
983 128/348 298 
983 476 - 
983 476/128/348 - 
462/521 128/348 ~ 150/150 
983/462/521 128/348 298/~150 


To test the reliability of the SNPs, amplicons of the five 
selected markers were digested with their corresponding 
restriction enzymes according to manufacturer's instruc- 
tions. We report here the conditions that we used, but these 
might need to be changed with the source of the enzymes, 
as different enzyme concentrations might require different 
digestion times. All digestions were performed at 37°C, and 
digestion solutions of TYRi, RAG1 and CMYC2 were incu- 
bated for 6 hours, while POMC and NTF3 digestions were 
incubated for 10 minutes. Digested fragments were then 
separated in 1% agarose electrophoresis gel for 30 min at 
100 V and 80 mA, and visualized with UV light. Fragment 
lengths were estimated by comparison with a 100-bp DNA 
Step Ladder (Promega). 

As explained above, we always aimed at genotyping 
each specimen for at least two diagnostic markers. When 
the RFLP patterns were ambiguous or discordant be- 
tween markers (each marker suggesting a different identi- 
fication), we repeated the digestion, and the amplification 
when necessary, until we obtained two unambiguous and 
concordant results. Samples that failed to produce unam- 
biguous RFLP patterns for at least two loci were discarded 
(almost invariably poorly preserved samples that failed to 
amplify reliably, see below). We only found two truly dis- 
cordant cases where clear RFLP results were consistently 
repeated on the same markers twice and produced RFLP 
patterns suggesting different species (see details below). 


Results and discussion 
Enzymatic restriction patterns 
on the markers selected 


In the results presented here, we report specimens whose 
identification was confirmed by previous allozyme geno- 
typing, by concordance between RFLP genotypings, or se- 
quencing of at least two diagnostic loci. Specimens that 
were genotyped for a single locus have been excluded and 
their identification has been treated as tentative only. 


RAGI1 with Drall 


The fragment amplified from the RAG1 marker (amplicon: 
983 bp, gene fragment of 936 bp after priming sequences are 


removed) contains several species-specific loci. Five diag- 
nostic SNPs differentiate P. ridibundus from P. perezi (two of 
which are shared between P. perezi and P. lessonae) and one 
SNP allele is specific of P. lessonae. We selected the Drall 
enzyme, with a target site of RGZGNCCY, that cuts between 
positions 565 and 566 for P. ridibundus alleles only (align- 
ment in Supplementary File 4), generating two fragments of 
589 bp and 394 bp for P. ridibundus but a single uncut frag- 
ment (983 bp) for P. perezi and P. lessonae (Table 3, Fig 1A). 
The digestions for P kl. grafi and P kl. esculentus produce 
three fragments, viz. one uncut band of 983 bp correspond- 
ing to the complete P. perezi or P. lessonae allele, and two 
smaller fragments of 589 bp and 394 bp, resulting from cut- 
ting the P ridibundus allele (Table 3, Fig. 1A, see Supplemen- 
tary File 5 for P lessonae and P. kl. esculentus patterns). 

This RAGi/Drall restriction site was genotyped by RFLP 
or sequencing in a total of 903 specimens (Supplementary 
File 1): 524 P. ridibundus (including three P. cf. ridibundus 
/ P. bedriagae from Azerbaijan), 86 P. perezi, 22 P. lessonae, 
251 P kl. grafi, and 20 P kl. esculentus; all yielded RFLP pat- 
terns or genotypes for the restriction site as was expected 
for their taxon (Table 3). 


TYRI with Ball 


The amplified TYR: fragment is a 646-bp amplicon (gen- 
erating a 601-bp fragment after removal of primer se- 
quences) where the enzyme Ball (target site: TGGzCCA) 
cuts between positions 318 and 319 in P. perezi, generating 
two fragments of almost the same length (305 and 341 bp), 
while in P ridibundus and P. lessonae this restriction site 
does not exist (alignment in Supplementary File 6). Asa re- 
sult, in the electrophoresis gel, a short fragment of around 
320 bp is visible for P. perezi, while a full-length fragment 
of 646 bp is found for P. ridibundus, P. lessonae and P. kl. es- 
culentus (Table 3, Fig 1B, see Supplementary File 5 for P les- 
sonae and P. kl. esculentus patterns). For P. kl. grafi, two 
bands of different length are visible, the shortest measuring 
around 320 bp, corresponding to the two fragments from 
the P perezi allele, and the longest 646 bp, corresponding 
to the complete P. ridibundus allele (Table 3, Fig 1B). Like 
in P. ridibundus, P. lessonae does not have a restriction site 
for the Ball enzyme, therefore digestion of P. lessonae and 
P. kl. esculentus samples results in a single band of 646 bp 
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(Table 3, see Supplementary File 5 for P. lessonae and P. kl. 
esculentus patterns). 

This TYRi/Ball restriction site was genotyped by RFLP 
or sequencing in a total of 886 specimens (Supplementary 
File 1): 523 P. ridibundus (including seven P. cf. ridibundus 
/ P. bedriagae from Azerbaijan and Turkey), 82 P. perezi, 13 
P lessonae, 251 P. kl. grafi and 16 P. kl. esculentus; all yielded 
RFLP patterns or genotypes for the restriction site as was 
expected for their taxon (Table 3). 


CMYC2 with SspI 


We did not manage to amplify the CMYC2 fragment in the 
PG system with the primers used. In the LE system, the 
amplicon is 298 bp, and the 253 bp fragment (after primers 
were removed) exhibits a unique species-specific SNP for 
P lessonae. The SspI enzyme, with a target site of AATZATT, 
only cuts the P lessonae allele, between positions 133 and 
134 (alignment in Supplementary File 7), producing a sin- 
gle 298-bp fragment for P. ridibundus and two undistin- 
guishable fragments of approx. 150 bp for P. lessonae (Ta- 
ble 3, Fig 1C). In P kl. esculentus, two bands of 298 bp and 
+ 150 bp are visible on gel (Table 3, Fig 1C). Even though 
we did not successfully amplify P. perezi, we retrieved the 
sequence of this species from GenBank, and as can be seen 


pag — 


eee 500bp 


300bp 


e 


R 


in the alignment in Supplementary File 7, P. perezi does 
not present the SspI restriction site. This implies that after 
digestion, the amplicons of P. perezi or P. kl. grafi would 
produce a pattern identical to P ridibundus and that the 
CMYC2z patterns described here would be truly diagnostic 
for P. lessonae and P. kl. esculentus. 

This CMYC2/SspI restriction site was genotyped by 
RFLP or sequencing in a total of 169 specimens (Supple- 
mentary File 1): 23 P. ridibundus, 78 P. lessonae and 67 P. kl. 
esculentus; all yielded RFLP patterns or genotypes for the 
restriction site as was expected for their taxon (Table 3). 


POMC with Hincll 


The fragment amplified from POMC (amplicon: 476 bp, 
415 bp of sequence after removal of primer) does not in- 
clude fixed differences between P. ridibundus and P. lesso- 
nae; however, it contains four SNPs that are only present in 
P. perezi. The HinclI enzyme, with target site GTYzRAC, 
digests the POMC fragment between positions 99 and 
100 (alignment in Supplementary File 8) in P. ridibundus 
and P. lessonae, generating fragments of 128 bp and 348 bp 
for these species and a full-length fragment (476 bp) for 
P perezi (Table 3, Fig 2). Among the hybridogenetic species, 
P kl. esculentus has the same restriction profile as its par- 


Figure 1. Enzymatic digestion patterns of the species diagnostic markers revealed by electrophoresis gel. Species names are labelled 
as follow. P: Pelophylax perezi; G: P. kl. grafi; R: P. ridibundus; L: P. lessonae; E: P. kl. esculentus. (A) Digestion of RAG1 by Dra II on 
species of the PG system (Digestion on the LE system is shown on Supplementary File 5). (B) Digestion of TYR1 by Ball on species 
of the PG system (Digestion on the LE system is shown on Supplementary File 5). (C) Digestion of CMYC2 by SspI on species of the 
LE system (Digestion on the PG system was not attempted). (D) Digestion of RAG1 by Bsm I on species of the LE system (Digestion 


on the PG system is shown on Supplementary File 5). 
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ents P. ridibundus and P. lessonae (two fragments of 128 and 
348 bp in length, Fig 2). Pelophylax kl. grafi presents three 
fragments, one of 476 bp, corresponding to the complete 
P perezi allele, and two smaller fragments of 128 and 348 bp, 
corresponding to the allele of P. ridibundus (Table 3, Fig 2). 

This POMC/HinclI restriction site was genotyped by 
RFLP or sequencing in a total of 83 specimens (Supple- 
mentary File 1): 18 P. ridibundus (including seven P. cf. 
ridibundus / P. bedriagae from Azerbaijan and Turkey), 19 
P perezi, 20 P. lessonae, 13 P. kl. grafi and 13 P. kl. esculentus; 
all yielded RFLP patterns or genotypes for the restriction 
site as was expected for their taxon (Table 3). 


RAGI1 with BsmI 


Initial examinations of RAG1 sequences suggested an ad- 
ditional diagnostic position at the enzyme BsmI restriction 
site CTTACzGN where this enzyme digests the P. lessonae 
RAGi allele between positions 437 and 438 (alignment in 
Supplementary File 4), generating two bands of similar 
length (462 and 521 bp) but usually visibly different (Ta- 
ble 3, Fig. 1D). For P. kl. esculentus, the digestion produces 
three fragments, one of 983 bp, corresponding to the com- 
plete P. ridibundus allele, and two bands of 462 and 521 bp 


Figure 2. Enzymatic digestion patterns of POMC by HmcII. DNA 
fragment size is determined by comparison with a 100-bp DNA 
Step Ladder (Promega). Species names are labeled as follow. 
P: Pelophylax perezi; G: P. kl. grafi; R: P. ridibundus; L: P. lesso- 
nae; E: P kl. esculentus. 


from the P lessonae allele (Table 3, Fig. 1D). For P. ridibun- 
dus, P. perezi and P. kl. grafi, the digestion results in a sin- 
gle band of 983 bp (Table 3, see Supplementary File 5 for 
P. perezi and P. kl. grafi patterns). 

As a consequence, a series of 146 specimens were geno- 
typed by RFLP for this RAG1/BsmlI restriction site (Supple- 
mentary File 1): 11 P. ridibundus, 1 P. perezi, 68 P. lessonae, 
1 P. kl. grafi and 64 P. kl. esculentus; all yielded RFLP pat- 
terns or genotypes for the restriction site as was expect- 
ed for their taxon (Table 3). However, when adding more 
sequence data, we realized that this restriction site is not 
diagnostic for P. lessonae bergeri from Italy, as two indi- 
viduals out of five sequenced Italian P lessonae carried the 
P. ridibundus allele at a heterozygous (BEV.11361) or ho- 
mozygous (BEV.12379) state (GenBank numbers in Sup- 
plementary File 1). These individuals would thus exhibit 
digestion patterns that do not correspond with the pat- 
tern expected for P. lessonae. Since we aimed at providing 
a method that would reliably identify the five widespread 
European species over their distribution range, we discard- 
ed the RAGi/BsmI for our final identification method, but 
this probably remains a valid option for molecular diagno- 
sis in most of Europe. Note however that DUFRESNEs et al. 
(2017a) reported an extensive occurrence of Italian alleles 
in most of the French populations of P. lessonae for two nu- 
clear loci so that the occurrence of Italian RAG alleles in 
French populations of P. lessonae, even if we did not detect 
them, is probably not unexpected. 


Repeatability of RFLP genotyping 


A variable but sometimes significant proportion of samples 
failed to amplify or yielded weak bands that were hard to re- 
solve on gel after digestion. This problem was mostly caused 
by poor PCR amplification due to low-quality DNA. Our 
samples originated from various sources (Supplementary 
File 1), including buccal swabs collected by non-profession- 
al naturalists and preserved either dried or in 95% ethanol, 
and old voucher specimens, generating samples of unequal 
quality. We especially noted that buccal swabs stored dry at 
room temperature were not as reliably amplified as those 
stored in ethanol. These issues are inherent to any PCR- 
based approach and not particular to our RFLP method. 

In addition, a small number of samples (45 out of 1,111, 
i.e., less than 5% of the genotyped samples) were success- 
fully amplified but yielded contradictory results between 
markers at a first attempt (genotyping results marked in 
yellow in Supplementary File 1). For most of these samples 
one or both markers were reamplified and/or the enzymat- 
ic digestion was repeated, resolving the contradictions in 
almost all cases. We suspect that most of these issues were 
caused by incomplete DNA digestion during the RFLP step, 
generating incorrect patterns. They were easily resolved by 
repeating the process (except for two samples, see below), 
however, this shows that the digestion step can sometimes 
generate flawed results that are impossible to detect if the 
identification is not confirmed by another marker. 
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Conflicting results between markers 


Two samples were consistently discordant between two 
markers even after repeated genotyping (orange results in 
Supplementary File 1): one sample from a mixed PG sys- 
tem and P. ridibundus population (T9844 in Supplemen- 
tary File 1) produced twice a P kl. grafi pattern for Drall 
on RAGI (371-, 565- and 936-bp bands) and a P ridibundus 
pattern for Ball on TYR1 (601-bp band); one sample from a 
LE population (T12991 in Supplementary File 1) produced 
twice a P kl. esculentus pattern for Bsm I on RAG1 (150- 
and 300-bp bands) and a P. ridibundus pattern for SspI on 
CMYC2 (300-bp band). The sample from the LE popula- 
tion presents P kl. esculentus morphological features and it 
was identified as such at sampling, but since morphological 
separation of P kl. grafi from P. perezi and P. ridibundus is 
difficult, we are unsure of the identification of the second 
individual. In these two cases the genetic identification was 
left unresolved, and we suspect that these contradictory re- 
sults stem from rare recombination events in the genomes 
of the hybridogenetic taxa (CHRISTIANSEN & REYER 2009). 
These rare cases of contradiction between markers, either 
due to recombination in hybrid taxa or allele-sharing at 
very low frequency in the parental species, are another rea- 
son that led us to design a method where identifications 
are always supported by a minimum of two independent 
markers. 


Re-examination of SAI-1 length polymorphism 


Havuswa pT et al. (2012) suggested that P lessonae differed 
from all other Western European Water Frogs by a ~530-bp 
deletion in SAI-1, allowing easy identification of P lessonae 
and P. kl. esculentus by the simple migration of the PCR 
products on agarose gel. We sequenced the SAI-1 fragment 
in two P. perezi individuals (GenBank OL702782, the other 
individual yielding a short partial sequence only). Together 
with the two specimens sequenced by Hauswa.pT et al. 
(2012, FN432377-78), these four sequences confirmed that 
P. perezi has a long SAI-1 allele, identical to P. ridibundus 
alleles after standard migration on agarose gels. As previ- 
ously suggested, SAI-1 length polymorphism should thus 
be an easy way of distinguishing P lessonae and P. kl. escu- 
lentus from the other Western European Water Frogs. 
However, when examining all (n = 19) P lessonae SAI-1 
sequences stored in GenBank, we found that three individ- 
uals (JQ965512, MFo94367-68) carry long SAI-1 haplotypes 
that cannot be separated from the alleles of the other spe- 
cies by migration on agarose gel. A quick examination of 
these sequences revealed that, while JQ965512 is identical 
to a widespread P. ridibundus allele and is possibly a mis- 
labelled sample, the other two sequences carry mutations 
that are private and group with other P. lessonae alleles in a 
neighbour-joining tree (result not shown). There is thus lit- 
tle doubt that P lessonae can carry long SAI-1 alleles, cast- 
ing doubt on the length polymorphism of this marker as 
a safe species-diagnostic tool in European water frogs. As 


226 


a consequence, we have not retained this marker for the 
identification method that we propose here. 


A novel PCR-RFLP protocol for the genetic identification 
of Western European water frogs 


We here propose a novel method to identify Western Eu- 
ropean Water Frog taxa by PCR-RFLP based on the re- 
striction patterns of selected markers. As we have detailed 
above, 1) the digestion step sometimes produces misleading 
RFLP patterns, and 2) there are very rare cases of discord- 
ance between identification suggested by different markers, 
therefore we developed a protocol that ensures that every 
identification result will be supported by a combination of 
two independent markers in all possible situations. 

Our method relies on three different combinations of 
markers, depending on the geographical situation. When 
only three species are present (P. ridibundus and species 
from either the LE system or the PG system), a combina- 
tion of PCR-RFLP tests on two markers will be enough to 
obtain a result. When all five species are potentially pre- 
sent, three markers are needed to ensure that all pairwise 
diagnoses are supported by two markers. Our protocol is 
detailed below and illustrated in Figure 3: 

1) in Spain and large areas of southern France, the only 
candidate species are P. ridibundus, P. perezi and P. kl. grafi, 
and a combination of RAG1/Drall & TYRı1/Ball facilitates 
reliable identification (Fig. 3A); 

2) in most of Europe, separation of P. ridibundus, P. lesso- 
nae and P. kl. esculentus can be achieved by a combination 
of RAG1/Drall with CMYC2/SspI (Fig. 3B); 

3) in the rare situations where all five species are poten- 
tially present (especially in western France), a first step with 
a combination of POMC/Hmcll with RAG1/Drall provides 
patterns that are unique to each species; a second step with 
TYR1/Ball for the PG system or CMYC2/SspI for the LE 
system allows to verify these identification results (Fig. 3C). 

With this protocol, all three pairwise identification re- 
sults between parental species are based on at least two 
markers: P. perezi and P. lessonae are separated by POMC 
and TYRu, P lessonae and P. ridibundus by RAGi and 
CMYCa2, and P perezi and P. ridibundus by RAG1, POMC 
and TYR1. We also made sure that no final identification 
result would be based on a combination of non-digested 
RFLP patterns to prevent misidentification due to com- 
plete digestion failure: when a species-specific RFLP pat- 
tern is “one band” (no digestion) for a given marker we 
always combine it with a marker that “cuts” for the same 
species. 


Applicability to other, exotic Water Frog taxa 


Although our protocol has been developed mainly on the 
basis of French samples, it has also been tested on a few 
samples from the LE system in Italy, where several parental 
and hybridogenetic taxa are currently treated as conspe- 
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cific with P lessonae and P. kl. esculentus (SPEYBROECK et 
al. 2020). When and where possible, we also checked on 
GenBank sequences that the substitutions used to design 
our tests were also present in Italian samples. Although in- 
creasing the sample size of Italian specimens would be ide- 
al, we believe that our method is also applicable to these 
populations. Within the P ridibundus complex, our data 
include specimens sampled in French populations carry- 
ing either the Central European mtDNA lineage or the Bal- 
kans (P. ridibundus kurtmuelleri) mtDNA lineage (based 
on data in DUFRESNES et al. 2017b and on our own data). 
As expected from their low divergence in nuclear DNA, 
these two lineages conform to the P. ridibundus patterns 
described above. We have not been able to assess the per- 
formance of our method in the P. ridibundus-bedriagae 
complex from Anatolia and the Middle Fast, but the sam- 
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ples that we have sequenced (Supplementary File 1, three 
from Azerbaijan and five from Turkey, belonging to several 
mtDNA lineages) share RAG1 and TYR: restriction sites 
with P. ridibundus from Europe, as has been assessed from 
our sequence alignments. 

Lastly, we wish to emphasize that our PCR-RFLP meth- 
od does not provide a “positive” identification result for 
any of the species: it distinguishes between these five spe- 
cies when the samples actually belong to one of these spe- 
cies. Samples of other species of the genus Pelophylax such 
as P. epeiroticus, P. sqhipericus, P. cretensis or P. saharicus 
will most likely amplify for the markers we used and could 
exhibit RFLP patterns identical to one of the five Western 
European species. For example, samples of the recently dis- 
covered introduced French population of P. saharicus (Do- 
NIOL-VALCROZE et al. 2021) were assayed for the RAG1/ 
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Figure 3. A novel PCR-RFLP protocol for the identification of Western European water frogs. (A) In the PG system, a combination 
of RAG 1/Drall and TYR1/Ball ensures that all identifications are supported by two independent loci. (B) In the LE system, a com- 
bination of RAG I/DralI and CMYC2/SspI has the same property. (C) When all five species are potentially present, the first step 
requires double digestion of POMC by HmclII and of RAG 1 by Drall, which in combination provides a species-specific pattern and 
reveals from which system (LE or PG) the samples are originating. This identification can be confirmed by a second step using TYR1/ 
Ball, if the first step suggests the PG system, or a CMYC2/SspI step, if the first step suggest LE. 
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Drall and TYR:/Ball combinations before their true iden- 
tity was revealed and they exhibited P perezi RFLP patterns 
for both markers (Supplementary File 9). Whenever our 
protocol returns an unexpected identification, either a spe- 
cies outside its normal distribution range (e.g., PG indi- 
viduals in Northern or Eastern Europe or LE individuals 
in Iberia) or an identification that does not match the mor- 
phology of the specimen, we urge that such identification 
be tested with other methods (e.g., DNA barcoding) to rule 
out the possibility of an introduced population of a species 
not included in our protocol. 


Conclusion 


Given the difficulties to identify European Water Frog spe- 
cies with morphological or acoustic characters and the con- 
servation concerns raised by the current spread of Pelophy- 
lax ridibundus at the expense of the native taxa of the PG 
and LE systems, DNA-based methods have become cru- 
cial for monitoring their status. The PCR-RFLP method we 
propose here generates restriction banding patterns specif- 
ic for P. ridibundus, P. perezi, P. lessonae and the two hybrid 
species P. kl. grafi and P. kl. esculentus, facilitating the reli- 
able identification of Water Frogs in most parts of Europe. 
In comparison to other molecular methods, our PCR-RFLP 
protocol is fast, relatively inexpensive, easily replicable, and 
allows non-invasive sampling. Because it systematically re- 
lies on at least two independent nuclear markers to reach 
an identification result, we also believe it is less error-prone 
than other PCR-based methods currently available. 
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